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Abstract 

Holographic Particle Image Velocimetry (HPIV) is a novel technique for measuring 
the complete fluid flow around a body. Advances in computing power make this 
technique practicable for the first time. 

Currently popular techniques for experimentally determining fluid flow around 
a test body rely on measuring the flow at a single point and moving the sample 
point during the experiment. This implicitly time averages the data. HPIV can 
measure the flow field in a volume nearly instantaneously by sequentially “holo¬ 
graphing” the seeded flow field. 

Until recently the only practical method of processing the HPIV data was to 
physically reconstruct the optical image from the hologram and scan through the 
resulting volume using microscope optics. This is a very slow process. 

Recent advances in computing power make it possible to directly analyze the 
digitized holograph (digitized using a high resolution scanner) This reduces the 
post-processing time from several hours to less than 1-hour with much less expen¬ 
sive equipment. 

The primary contribution of this thesis is an advanced multi-platform graphical 
user interface hologram synthesizer for use in calibrating digital HPIV reconstruc¬ 
tion algorithms. 

Thesis Supervisor: Jerome H. Milgram 
Title: Professor of Ocean Engineering 
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Chapter 1 
Introduction 

1.1 Overview 

This thesis presents an advanced multiplatform graphical user interface hologram 
synthesizer to produce calibration holograms of seeded flow fields. The synthe¬ 
sizer, HoloSynth, runs under X11R6 with OpenGL, and Microsoft Windows NT. 
Complete source code is contained in Appendix A. A brief introduction to parti¬ 
cle holography and holographic particle velocimetry is presented as well as results 
from computational reconstruction of the hologram generated by the synthesizer. 

1.2 Motivation 

Until very recently, particle image velocimetry (PIV) was limited to measuring 
fluid flow velocities only in a single plane. A sheet of laser light was passed 
through the flow which contained neutrally buoyant particles. A double exposure 
of the particles onto either film or a Charge Coupled Device (CCD) camera allows 
the two velocity components in the illuminated plane to be determined. 

The PIV method has recently been extended to three dimensions by Katz et. al. 
[1] by making a hologram of the flow field. Three dimensional Holographic Particle 
Image Velocimetry (HPIV) uses either in-line or off-axis holographic techniques 
(see Section 2.2). Once the hologram has been recorded, the image is reconstructed 
by reversing the recording process and using a three dimensional optical scanner to 
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digitize the image for analysis by computer. This process can take several hours. 

An alternative to scanning the reconstructed image is to digitize the hologram 
and allow a computer to reconstruct the image in software. The recent advances in 
computing power should make this process far faster than the mechanical/optical 
scanning process now used. 

Digitizing the hologram can be done in two ways. The initial method avoids 
the three dimensional scanning of the reconstructed image by using a high perfor¬ 
mance flat bed scanner to digitize an enlargement of the transparency produced 
in the normal holographic process. Alternately, a high resolution scanner can be 
used directly on the hologram transparency. The flat bed scanning process is well 
understood, since the transparency can be enlarged to overcome resolution prob¬ 
lems with the scanner. However, it still requires recording the hologram on film 
with the attendant film processing. 

A more advanced technique would be to capture the hologram using a Charge 
Coupled Device (CCD) Camera. The hologram would be immediately available 
in digital form, and no film or processing would be involved. 

Regardless of the hologram digitization technique the computational machinery 
for reconstructing the image must be calibrated. It far easier to produce a synthetic 
hologram for calibration than to set up a tightly controlled experiment and actually 
record a hologram. 


1.3 Requirements 

The initial goals of computationally reconstructed HPIV are to be able to accu¬ 
rately determine the flow field in a cubic volume a few inches on a side. The fluid 
speeds will be typical of those used in hydrodynamics research. The quantifiable 
requirements are; 

• Must be able to resolve a 10/um particle. This corresponds to the most 
conveniently sized neutrally buoyant seed particles used for Laser Doppler 
Anemometry. 
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• Accurate coverage within a 125cm 3 (5cm cube) volume. This will rise as the 
technique is perfected, and computing power increases. 

• Unambiguously detect fluid speeds < lOm/s. This accounts for free stream 
plus allowance for turbulence and model scale factors. 

The first and second requirements dictate the size and resolution of the recording 
device. 

Since holography is essentially a linear process, multiple exposures can be used 
instead of a rapid succession of exposures. So rather than design a system to 
rapidly capture two separate holograms, the recording medium is exposed two 
or more times without being altered between exposures. This results in multiple 
real images for each particle in the volume, the analysis process must identify the 
multiple images with a particular particle. This alleviates the need for a high 
speed film traversal mechanism and removes CCD output data rate as a concern. 

The third requirement sets the pulse lengths and pulse intervals of the laser 
used to generate the hologram. If the pulses are too long the diffraction patterns 
will be smeared from the particle motion. The greater the interval between pulses, 
the harder it will be to uniquely identify particle images with a single particle. The 
pulse length and interval requirements determine the recording device sensitivity 
and laser output power. 
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Chapter 2 

Review of Holographic Technique 

2.1 Overview 

A hologram is a two dimensional recording that can produce a three dimensional 
image. Typically electromagnetic (light) waves are recorded, but recording of 
acoustic waves and matter waves have also been made[2]. Electromagnetic holog¬ 
raphy differs from normal photography in that it encodes phase information as 
well as intensity information on the recording medium. 

Recording media respond only to the intensity of the incoming wavefront, not 
its phase. To encode the phase using normal media, holography depends on using a 
coherent reference beam to interfere with a scattered or diffracted wave, or “object 
wave”, from the objects under observation. As a result of this interference, the 
intensity at the recording plane depends on both the phase and amplitude of the 
object wave. 

The original image is reproduced simply by illuminating the hologram with 
coherent light using geometry similar to that used to record the hologram. If the 
wavelength used for reconstruction differs from the wavelength used for recording, 
the image will be enlarged or reduced. 

The reconstruction process, while simple, is not amenable to rapid quantitative 
data collection. Quantitative analysis generally involves using three dimensional 
optical scanning of the reconstructed image. The ability to computationally re- 
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Figure 2-1: Basic in-line hologram geometry 


construct the image would dramatically increase the utility of holographic mea¬ 
surement techniques. 


2.2 In-line (Gabor) Holography 

In-line holography, or Gabor Holography after its inventor, describes a holographic 
apparatus in which the reference beam is usually oriented perpendicular to the 
recording media and illuminates the objects from behind. It is not a useful tech¬ 
nique for objects that are large with respect to the recording medium, but if one 
is interested primarily in the three dimensional position of small objects and not 
their three dimensional shape, in-line holography can be very effective. It is com¬ 
monly used to size particulates and aerosols suspended in fluids, and it has been 
used for traditionally reconstructed HPIV[3]. 

Figure 2-1 shows the basic geometry and coordinate system for recording a 
plane in-line hologram. The object is illuminated by the reference source. The 
diffracted wave, or object wave, interferes with the wave from the reference source, 
or reference wave. This interference pattern is recorded by the medium in the 
recording plane and constitutes the hologram itself. 

If the reference and object complex wave amplitudes measured at the recording 
plane are given by R ( x , y) and O (x, y) respectively, the intensity at the recording 
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plane is given by, 


I(x,y) = \R(x,y) + 0(x,y)\ 2 (2.1) 

Since the recording plane is perpendicular to the plane reference wave, the 
complex amplitude of the reference wave is a constant, r. Expanding (2.1) yields 

I(x,y) = |r + O (x, y)\ 2 

= r 2 + \0(x,y)\ 2 + rO(x,y) + rO* (x,y) (2.2) 


The intensity is recorded in any acceptable manner, usually as a positive trans¬ 
parency on film. Most reasonable recording media respond nearly linearly to in¬ 
cident wave intensity. The transmittance, t(x,y), of a positive transparency can 
be adequately described by 


t(x,y) = r b — pTI (x,y) 

= n - f3T jr 2 + |O (x , y)\ 2 + rO (x, y) + rO* (x, y) 


(2.3) 


where r b is a constant background transmittance, T is the exposure time, and /3 
is a constant of the recording medium. 

Figure 2-2 shows an in-line hologram being physically reconstructed. The 
wavefront transmitted by the hologram is simply the incident wave multiplied by 
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the transmittance of the hologram and can be written, 

ip(x,y) = rt(x,y) 

- r (r b +/3Tr 2 ) + PTr\0(x,y)\ 2 

+/3Tr 2 0 (.x , y) + (3Tr 2 0* (x, y) (2.4) 

The first term in (2.4) represents the uniformly attenuated reference wave. The 
second term is normally very small, since typically O ( x , y) <C r. It is called the 
“halo” and under proper conditions can be neglected. The third term is identical 
to, within a multiplicative constant, with the original scattered wave. It forms an 
image z 0 from the hologram on the side opposite the observer. This is the virtual 
image, which is what the observer perceives. The fourth term represents a wave- 
front similar to the originally scattered wavefront, but with opposite curvature. It 
converges to form a real image z 0 from the hologram, between the hologram and 
the observer. [4] 

The observer then sees the virtual image, with a superimposed, out-of-focus 
real image along with the diminished reference beam. The two images instead of 
one and are the biggest drawbacks to in-line holography. 


2.3 Off-axis (Leith-Upatneiks) Holography 

A simple way to separate the out-of-focus real image from the virtual image is to 
provide a separate path for the reference beam and the object wave. Figure 2-3 
shows the basic geometry for recording an off-axis hologram. 

Because the reference beam is no longer perpendicular to the recording plane, 
it cannot be simplified to a constant, rather, it can be expressed as, 

R (x, y) = r exp(j2n£, r x) (2.5) 

where = (sin#)/A, since only the phase of the reference beam varies across the 

recording plane. 
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In a fashion similar to the derivation of (2.2) the intensity at the recording 
plane is found to be, 

I(x,y) = r 2 \0(x,y)\ 2 + 2r |0 (x,^)| cos [27r£ r x + </> (x,y)] (2.6) 

Figure 2-4 shows the geometry for an off-axis hologram reconstruction. As¬ 
suming the recording medium is linear, the transmitted wave complex amplitude 
is found in a similar manner, 

ip(x,y) = r b r exp(j2n£ r x) 

+PTr |O (x, y) | 2 exp(j27r£ r x) 

+/3Tr 2 0(x,y ) 

+(3Tr 2 0 (x, y) exp(j47r£ r x) (2.7) 

Once again, the first term in (2.7) is simply the attenuated reference beam. 
The second term is a “halo” term, similar to the neglected term in (2.7). The 
third term is identical to the original scattered wave and forms a virtual image 
in the same position as the original object. The fourth term is, again, the real 
image, this time shifted from the optical axis by approximately twice the angle 
the reference beam makes with the recording plane. Hence the out of focus real 
image does not interfere with the vitual image. 
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At first blush it would seem off-axis holography is the answer to the virtual 
image problem for HPIV. If purely physical reconstruction is to be used, then off- 
axis holography is eminently practicable. The fatal flaw comes when considering 
digitizing the hologram for use in computational reconstruction. The reference 
beam now introduces a rapidly fluctuating intensity term that quadruples the 
resolution requirement of the recording medium and digitization system. For this 
reason off-axis holography will not be discussed further in this thesis. 


2.4 Diffraction 

Holography is fundamentally about diffraction. Goodman[5] provides a detailed 
derivation of the Huygens-Fresnel principle for a wave from constant 

Tp(xh,yh) = JJ ipz{x, y) eXP ^ r °^ dx dy . (2.8) 

£ 

where, r 0 i = \/z 2 + (xh — x) 2 + — y) 2 , and Xh and y h are coordinates in the 

recording plane. 

This expression, while physically meaningful, is not easily evaluated. Two 
fundamental simplifications are common depending on the geometry of the optical 
problem being analyzed. The approximations are the Fresnel and the Fraunhofer 
approximations, valid in the near and far fields, respectively. 
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2.4.1 The Fresnel Approximation 

A simpler and more useful expression of the Huygens-Fresnel integral can be found 
by simplifying r 0 i using the paraxial approximation, z 2 x 2 + y 2 . 

If b < 1, then, 

\/l + b = 1 + — b — —b 2 + • • • (2-9) 

2 o 

Using (2.9) with b = ( ~ x ) 2 + 


r oi = \[z 2 + (x h - x ) 2 + ( y h - y ) 2 


= zWl + 


x h - X 


+ 


j/fc-y 




( 2 , 10 ) 


Substitution of (2.10) into (2.8) yields little simplification unless carefully cho¬ 
sen terms are dropped. For the term in the denominator of (2.8) it is considered 
acceptable to drop all terms but 2 . However, small changes in the roi occurring 
in the numerator can cause significant changes in the value of the exponential 
and cannot be omitted. These large changes are generated in part by the mul¬ 
tiplication of r 0 1 by k, usually a very large quantity, and the sensitivity of the 
exponential to phase changes in its argument. The resulting simplification yields 
the Fresnel-Kirchoff Diffraction Integral: 

jkz °° 

^(x h ,y h ) = JJ U(x,y)exp^j^ [(x h -x ) 2 + (y h - yf]^j dx d&.ll) 

— OO 

jkz r 

= Jj U(x,y)e j ^( x2+y2 )e- j ^ XhX+ v^dxdy( 2 . 12 ) 

—OO 


Where (2.12) differs from (2.11) only through factoring the term exp {x\ + y 2 )] 
outside the integral. This approximation holds when the observer is in the near 
field, or Fresnel zone, of the diffracting aperture. Typical HPIV geometry will 
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place the recording medium in the far field of the individual particles but in the 
near field of the ensemble. 


2.4.2 The Fraunhofer Approximation 

Additional simplifications to the Fresnel Diffraction Integral can be made if the 
system meets the Fraunhofer approximation, 


k(x d - y )max 
2 > -~- 


(2.13) 


Simplifying (2.12) yields, 


il’(xh,yh) = JJ U(x,y)e i^ XhX+yhy ' 1 dxdy (2.14) 


which is valid in the far field. The integral is recognizable as a two dimensional 
Fourier transform. 


2.4.3 Convolution interpretation of Fresnel-Kirchoff Inte¬ 
gral 

The Fresnel-Kirchoff Integral, 


ip(x h 


fVh) = JJ U(x,y)ex p j [( x h - xf + (y h - yf] J dxdy (2.15) 


can easily be viewed as a convolution with kernel, 


oikz 


h z (x,y) = — exp 




(2,16) 


yielding 


V’s 


OO 

(Xh,yh) = JJ u(x,y)h z (x h - x,y h - y)dxdy = U(x h ,y h ) ® h z (x,y) (2.17) 
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This expression determines the field diffracted from a planar object, or an ensemble 
of co-planar planar objects. If a large number of very small non-coplanar objects 
are in the field of view, (2.17) can be extended: 

1>(xh,Vh) = ^2u n (x h ,y h ) ®h Zn (x h ,y h ) (2.18) 

n 

where the sum is taken over n planes within the volume. Each object plane, 
Unix, y), contains at least one diffracting object and is located z n from the record¬ 
ing plane. The planes need not be evenly spaced. 

2.5 Laser Characteristics 

2.5.1 Fundamentals 

At the atomic level, light is always generated by an atom transitioning from one 
energy state, E 0 , to a lower energy state, E \. When this transition occurs the atom 
emits a photon with frequency v = (E 0 — Ei)/h, where h is Planck’s constant. 
This emission is a probabilistic occurrence and may occur spontaneously provided 
the atom has sufficient internal energy. This is spontaneous emission. 

A second form of emission, first suggested by Einstein in 1917, is called stim¬ 
ulated emission. In the case of stimulated emission, a photon collides with an 
excited atom. The incident photon must have exactly the same energy as the 
transition the atom will undergo. When this condition is met, the atom will make 
the transition, emitting a second photon. The emitted photon is in phase with 
the incident photon and travels in the same direction. 

Stimulated emission is fundamental to laser operation. A material is induced to 
have a majority of its atoms in an excited state (called a “population inversion”). 
Since the majority of the atoms in the material are excited to the same state, 
they have identical transition energies. Initial spontaneous emission by a small 
number of the atoms will stimulate emission in others. The stimulated photons 
will in turn stimulate further emission in a run-away chain reaction. This in itself 
will not cause light amplification, but if the material is mounted within an optical 
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cavity (two parallel reflectors), the stimulated emission will preferentially orient 
itself along the axis of the cavity. Radiation not parallel to the axis will be lost 
to the surroundings and not contribute to the electromagnetic wave building up 
in the cavity. Since the stimulated emission is always in phase, the radiation 
traveling along the cavity axis constructively reinforces and the amplitude rapidly 
increases forming a longitudinal vibration mode within the cavity. If an end of 
the cavity is only partially reflective, then some of the coherent amplified signal 
will be released, this is the usable laser beam. 

2.5.2 Establishing the Population Inversion 

The process of establishing the population inversion is called ‘pumping’. Pumping 
can be accomplished using flash lamps, usual in the case of ruby lasers, or electrical 
discharge, as in the case of most gas lasers. Photochemical reactions can also 
stimulate the population inversion but that technique is very rare. [6] 

The cavity material characteristics fundamentally determine what type of pump¬ 
ing will be appropriate. Continuous beam lasers must obviously use a sustainable 
pumping technique, such as tungsten filament incandescent lamps, or the afore¬ 
mentioned electrical discharge. 

The length of the pulse in non-switched pulse lasers is directly related to the 
length of the pulse length of the pump. Flash lamp pumped ruby lasers typically 
provide pulses on the order of several milliseconds long. The length of the pulse 
also varies from pulse to pulse as the condition of the flash lamp changes. 

The minimum time between lasers pulses in non-switched lasers is directly 
related to the requirements of the pump. Flash lamps typically must cool imme¬ 
diately after use. The amount of cooling is directly related to the power output 
of the lamp, and can be on the order of minutes. This limits the pulse repetition 
frequency of the unswitched laser. 

Switching can reduce pulse lengths to femtoseconds (lxlO -15 seconds) and 
pulse intervals similarly. Switching is discussed in some detail in Section 2.5.5. 
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2.5.3 Energy Levels 

When common laser materials are pumped, the resulting atomic population usu¬ 
ally has many more than one possible transition state available. The likelihood 
that an atom undergoes a particular transition depends on the quantum mechan¬ 
ics of the particular material. In any event, when the population inversion starts 
to break down there can be a large number of longitudinal oscillation modes, 
each producing a separate wavelength in the cavity. Typically, a small number of 
closely spaced wavelengths will predominate, but the resultant beam is not strictly 
monochromatic without special modifications to the laser (see Section 2.5.7). [6] 

2.5.4 Coherence Length 

The preceding discussion leaves the impression that all light generated in the laser 
cavity is coherent. This is not the case. Not all stimulated photons are descended 
from the same spontaneously emitted photon, and the Heisenberg uncertainty 
principle places a lower limit on the precision with which the wavelengths are 
reproduced. So even a laser oscillating in only a single longitudinal mode has 
some degree of incoherence in its output. 

One measure of the quality of a laser’s output is termed “coherence length”. 
This length is the maximum beam path length within which the beam can coher¬ 
ently interfere with itself. The practical ramification of this is that the holographic 
system must have all path lengths from laser output to recording planes less than 
the coherence length of the laser. [7] 

2.5.5 Switching 

It is clear that the amplification properties of the laser cavity are critically depen¬ 
dent on the alignment and quality of the end mirrors. Misalignment or surface 
imperfections can destroy the amplification effect and turn the laser into an ex¬ 
pensive flash bulb. 

This sensitivity can be turned to the laser users advantage. If the mirror prop¬ 
erties can be accurately controlled then the amplification process can be controlled 
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independently of the pumping process. 

Borrowing a term from electrical engineering, the ability of the cavity to res¬ 
onate is called it’s ‘Q\ Systems that switch the amplification process on and off 
by controlling the resonance qualities of the cavity are called ‘Q-switched’ or ‘Q- 
spoiled’ lasers. 

Q-switching has a subtle benefit. The population inversion does not occur 
instantaneously at the beginning of the pumping sequence. If the cavity cannot 
oscillate due to low Q, the population inversion can build to a much higher level 
than required for lasing. If the Q is suddenly increased, allowing the cavity to 
oscillate, the peak energy in the pulse can be much higher than if the cavity was 
allowed to oscillate continuously through the pumping sequence. So the power 
available in a Q-switched pulse can be much higher than in an unswitched pulse. [6] 

An additional benefit of Q-switching is the ability to produce multiple pulses 
out of a single pumping process. This allows accurate control of continuously 
pumped lasers, and multiple rapid pulses out of short pulse pumped lasers such 
as ruby lasers. 

2.5.5.1 Electrooptical Q-Switching 

Several physical phenomena allow the electrical modulation of an optical signal. 
The most commonly used in Q-switching lasers is embodied in the Pockels Cell. 

Pockels Cells are formed from a crystal, typically potassium dihydrogen phos¬ 
phate (KDP). KDP has two distinct optical axes. With an electrical voltage 
applied across the crystal the index of refraction on one axis is different than the 
other. An incident polarized wave is split into two portions, one portion traveling 
more slowly than the other. The two waves are recombined upon exiting the cell, 
resulting in a wave that is elliptically, circularly, or linearly polarized depending 
on the applied voltage. 

The two voltages of interest for Pockels Cell Q-switching are termed the quarter- 
wavelength (A/4) voltage, and the half-wavelength (A/2) voltage. The quarter- 
wavelength voltage will result in a circularly polarized beam at the output of the 
cell. The half-wavelength voltage will result in a wave linearly polarized 90° to 



2.5. LASER CHARACTERISTICS 


24 


Output Mirror 
Laser Rod 



Output Mirror 


Laser Rod 


Polarizer 

Pockels Cells 

Rear Mirror 



A/4 
Retardation 
induced by 
applied voltage 



(a) 



Polarizer 1 


Pockels Cells 


A/2 

Retardation 
induced by 
applied voltage 


_^^Polarizer 2 



Rear Mirror 

(b) 


Figure 2-5: Pockels Cell Q-switch at (a) quarter-wave and (b) half wave retardation 
voltage. 


the input wave. 

The quarter-wavelength retardation system shown in Figure 2-5(a) is a ‘pulse 
off’ Q-switch. During the initial phase of pumping, the A/4 voltage is applied to 
the cell. The output of the Pockels cell is circularly polarized, after reflection the 
circularly polarized beam reenters the Pockels cell and emerges linearly polarized 
90° to the input beam. The polarizer then prevents the beam from reentering 
the laser rod, spoiling the amplification. When the voltage is removed, the Q is 
restored, and oscillation can start. 

The half-wavelength retardation system shown in Figure 2-5(b) is a ‘pulse on’ 
Q-switch. With no voltage applied the Pockels cell does not alter the polariza- 
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Figure 2-6: Schematic Cavity Dump 


tion of the transmitted beam, and the crossed polarizers prevent oscillation from 
occurring. When switched on the Pockels cell shifts the polarization 90° allowing 
the beam to pass uninterrupted through the two polarizers. 

The half-wave retardation voltage for KDP is 7500V. KDP is also water soluble. 
Normally the crystal is hermetically sealed to prevent exposure to atmospheric 
moisture, but hydrophobic materials and extremely high voltages merit caution 
when used in a hydrodynamics experiment. 

2.5.5.2 Mechanical Q-Switches 

There are two common methods of mechanically Q-switching a laser. They both 
rely on mechanically rotating a reflective element within the laser cavity. The most 
obvious is to rotate the mirror such that it is only periodically parallel to the fixed 
reflector at the opposite end of the cavity. This system requires extremely critical 
alignment of the shaft upon which the mirror rotates. The alternative to rotating 
a mirror is to use a roof prism, which alleviates the alignment requirement. 

Mechanical switching is not as flexible as electrooptical switching, and is prob¬ 
ably not practicably capable of implementing the three pulse technique described 
in Section 3.1. 
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Figure 2-7: Cylindrical Transverse Electromagnetic Modes 


2.5.5.3 Cavity Dumping 

Extremely short pulses can be generated using a “cavity dump”. Cavity dumps 
allows the use of two 100% reflectivity end mirrors. When the laser power is 
maximized, the cavity dump redirects the beam outside of the cavity, rapidly 
dumping the entire optical energy of the cavity into the output beam. The pulse 
length is reduced to the round trip transit time within the cavity. Cavity dumps 
can be used to switch intermittently and continuously pumped lasers. 

A schematic of a cavity dumped ruby laser is shown in Figure 2-6. Initially, no 
voltage is applied to the Pockels cell. When the flash lamp is fired the horizontally 
polarized wave is passed through the polarizer and not allowed to regenerate. 
Vertically polarized light is kept within the cavity and builds up the vibration mode 
within the cavity. When the energy stored in the Ruby has peaked, the Pockels 
cell is switched to its half-wave retardation voltage. The vertically polarized light 
is rotated to horizontally polarized light by the Pockels cell and passes through the 
polarizer and out of the cavity. The combination of the polarizer, off-axis mirror 
and the Pockels cell essentially forms a voltage controlled mirror, whose reflectivity 
can be change from 100% to 0% and back to 100% during the cycle. [6, 8] 
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2.5.6 Transverse Electromagnetic Modes 

Laser cavities do not only oscillate longitudinally. There are also transverse modes 
of vibration which are very similar to the vibration modes of a drum head. Laser 
specifications typically include a description of which transverse modes are present 
in the output. The Transverse Electromagnetic Modes (TEM) are described using 
the symbol TEM mng . The first two indices describe the number of transverse 
modes, while the third (often omitted in specification) describes the number of 
longitudinal modes. Figure 2-7 shows cylindrical transverse modes, along with the 
indices describing them. Rectangular modes are also possible. It is not difficult 
to find TEMoo ruby lasers suitable for HPIV use. 

2.5.7 Output Frequency Control 

It is possible to dramatically reduce the number of longitudinal modes present in 
the output. The simplest method is to use a Fabry-Perot etalon. This is simply 
a very small optical cavity that resonates at precisely the frequency desired in 
the laser cavity. The etalon effectively spoils the cavity Q for other frequencies, 
preventing them.from building up. 

The simplest form of etalon is simply a thin glass plate. Internal reflections 
within the plate selection interfere with the wave to produce a sharp resonance 
peak. The effective thickness of the plate can be changed by rotating the etalon 
about an axis perpendicular to the cavity axis. 

Other etalon designs are possible. Alternate designs include wedge shaped 
plates that can be tuned by translating the plate, as well as designs that form 
an optical cavity using an air gap between two plates. The latter design can 
sometimes be adjusted by varying the gas pressure or temperature within the 

J 

cavity. [6, 8] 

2.5.8 Laser Construction 

Commercial lasers for field use are often manufactured in such a way as to limit the 
ability of the user to alter the laser. This increases laser reliability. Commercial 
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lasers for laboratory use are a different matter. As discussed above there are 
many additional components that can be added to the laser cavity to customize 
and control the laser output. Laboratory lasers typically mount all components, 
including the mirrors and laser rod, on a long rail with additional space for adding 
polarizers, etalons, Q-switches and cavity dumps. 
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Chapter 3 

HPIV System Design 


3.1 Timing 

The requirements to unambiguously determine speeds less than lOm/s inside a 
volume 5 cm on a side impose severe restrictions: 

1. Exposure times must be short enough to ‘stop’ the motion. In a photograph 
the motion allowed during the exposure is related to the feature size of the 
object. During the photographic exposure the smallest features of an object 
should not move more than approximately 1 /10th their size or the image will 
be unacceptably blurred. Likewise an acceptable amount of particle motion 
during holographic exposure has been found to be approximately l/10th of 
a particle diameter [9]. 

2. The interval between exposures must be short enough to maintain the iden¬ 
tity of the particles, so tracking and velocity determination is possible. This 
restricts the mean free path length of a particle between exposures to less 
than the distance between particles. There is an additional problem of de¬ 
termining which particle image is ‘oldest’, so that direction of travel can be 
determined. 

3. The interval between exposures must be short enough to ensure a majority 
of the particles captured in the first exposure are within the covered volume 
during the following exposures. 
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For a 10/um particle, these requirements become: 

1. A 1/us exposure results in l/10th diameter of travel during the exposure. 

2. A 10/us exposure interval results in 10 particle diameters of travel between 
exposures at lOm/s. 

The exposure interval also places an upper limit on the concentration of seed 
particles. Assuring the particle mean free path between exposures is less than the 
average distance between particles requires that 


r\ 


Where p p is the seed particle density, and the mean free path of the particles is 
ten particle diameters, 100/um. 

The direction problem is subtler. Keeping the average inter-particle distance 
greater than the mean free path length between exposures significantly simplifies 
the task of identifying pairs of particle images as being the same particle. Unfortu¬ 
nately, even when the two images of a given particle have been identified, there is a 
180° ambiguity in the velocity of the particle. A simple solution to this problem is 
to use three pulses rather than two, with differing time intervals between pairs of 
pulses. For example, the second 1/us pulse is sent 5/us after the first, then the third 
is sent 10/us after the second. The pulse intervals are sufficiently short to permit 
particle identification, while the difference in intervals allows determination of the 
time order of the images. An additional advantage of the triple exposure tech¬ 
nique could be estimation of fluid acceleration. This technique assumes the flow 
speed will not change appreciably over the distances traveled by particle during 
the exposure sequence. If the particle were to radically decelerate to less than half 
of the speed it traveled during the first set of exposures, it would be possible to 
confuse the time order of the exposures. 
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3.2 Recording and Digitization 

3.2.1 CCD Requirements 

The ultimate goal for computationally reconstructed HPIV is a completely digital 
process for determining the flow field throughout the volume. The most likely can¬ 
didate for the ‘digital recording plane’ is a Charge Coupled Device (CCD) Array. 
Unfortunately, the resolution of CCD arrays is significantly less than holographic 
film. 

3.2.1.1 Characteristics of CCD arrays 

A CCD Array is a semiconductor architecture that transfers charge through stor¬ 
age areas. The CCD array is essentially an array of charge “buckets”. Photons 
are captured by the buckets and converted to charge, the buckets only perform 
this conversion during a controllable integration period. More photons, or higher 
energy photons, hitting the array during the integration period result in a higher 
element charge. After integration the CCD Array support electronics scan through 
each bucket and report the charge level in each bucket. This charge is proportional 
to the intensity of the image focused on the array. [10] 

For holographic applications the fundamental characteristic of CCD Arrays 
that differentiates them from , film is the limited number of recording elements. 
High quality photographic film can have grains smaller than 1/rm over an area 12- 
20 centimeters square, or better than 120,000x120,000 recording elements. The 
best CCD arrays currently commercially available have element sizes on the order 
of 10/im, but the total array size is only about 5 centimeters on a side, yielding 
an array of approximately 5000x5000 elements. This is the root of the limitation 
in holographic applications. 

3.2.1.2 CCD Resolution Requirements 

3.2.1.2.1 Diffraction patterns of circular objects Parrent and Thomp¬ 
son analytically derived the diffraction pattern for an opaque circular object in a 
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coherent field. [11] For a circular particle, the field and intensity, respectively, are: 



Where, r is measured from the center of the object projected into the recording 
plane, and A(a) = Ji(q)/q, a Bessel function of the first kind divided by its 
argument. 

An illustrative example is plotted in Figures 3-1 and 3-2. The parameters 
more characteristic of HPIV systems result in plots that do not acceptably show 
the salient features of the intensity distribution. 


3.2.1.2.2 Recording the Diffraction Pattern Equation (3.1) describes the 
diffraction pattern recorded for a single particle with circular cross section. All 
previous discussion has assumed this diffraction pattern is perfectly recorded. In 
reality the recording medium has limited bandwidth ( i. e. a grain), and limited 
physical extent. 

Figure 3-1 shows a high frequency sine function modulated by a lower frequency 
A(a) envelope function. The envelope function is dependent on both particle 
diameter and particle distance from the recording plane. The high frequency 
sine function is dependent only on particle distance. Accurately recording the 
high frequency component is the most critical factor in obtaining good depth 
discrimination. 

Vikram[3] recommends recording the first three side lobes formed by the bessel 
envelope function. The spacing of the fine fringes, Ax, in the third lobe can be 
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r (cm) 

Figure 3-2: Enlarged plot of the third side lobe in Figure 3-1 


found using the argument of the sine term in (3.1), kr 2 /2z: 

u„ 2 

n/dyl 

--= Z7T 

2 z 2 z 

x\ — x\ ~ 2 xAx = 2 Az 

Ax = (3.2) 

x 

Equation (3.2) should be evaluated in the region containing the third side lobe. 
The mid-radius of the third side lobe, f 3 is: 


Ji(a n )/a n = 0 (n = 1 , 2 ,..., oo) 


&n. — 


T n = 


rs = 


a n = 3.832, 7.06, 10.173, 13.323, 16.470. 
kar n 
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Tn T T n —i 


2ka 


(o^Ti + ®n—l) 


( 0:3 + 0 : 2 ) = 17.189 


2ka 


2 ka 


(3,3) 


Equations (3.3) and ( 3 . 2 ) give f 3 = 0.189cm and Ax = 3.65 x 10 -3 cm for 
the parameters plotted in Figures 3-1 and 3 - 2 . Using z = 5 cm and a = 10 /mi, 
with A = 690nm, more reasonable for HPIV purposes, r 3 = 0.472cm and Ax = 
7.31 x 10 _ 4 cm. 

The parameter f 3 dictates the minimum size of the hologram, with no mag¬ 
nification, and is easily met. The parameter Ax dictates the minimum grain of 
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the hologram recording medium. Since enlargement or reduction of the hologram 
during recording is easily accomplished using lenses, the absolute size of the fringe 
spacing is not important. Rather, the important quantity is the ratio the fringe 
spacing to the total size of the system, so that the fine fringes produced by par¬ 
ticles at opposite extremes of the volume can be resolved. The sampling theorem 
suggests that the recording medium should be able to sample at twice the high¬ 
est spatial frequency required, or 1/2 the fringe spacing. To achieve the desired 
accuracy, a CCD array would have to have 

2 • 0.05cm . , 

——-—— : — « 13680 pixels square. 

7.31 x 10 _4 cm 

This resolution will adequately record the particle diffraction patterns. Unfortu¬ 
nately no CCD arrays are commercially available that reach this resolution. The 
best currently available are in the range of 5000 pixels square. 

3.2.2 Exposure Requirements 

While the performance characteristics of CCD arrays are well known it is difficult 
to predict the amount of power incident on the array from the suspended particles. 

The resolution analysis was performed assuming the finest necessary fringes 
were detected. The incident power depends heavily on the power output of the 
laser, the absorption of the fluid, and the reflectivity of the particles, and the 
viewing window material. 

The exposure problem is different in character from the resolution problem. 
The resolution problem has to wait for better CCD technology, or lower the ac¬ 
curacy requirements for HPIV significantly. The exposure problem can be solved 
many ways, the simplest (but not necessarily the least expensive) way is to increase 
the laser power. 
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Chapter 4 
Simulation 

4.1 Motivation 

A primary factor preventing HPIV from becoming an easily and widely used tech¬ 
nique is the data reduction effort. As previously mentioned, physical reconstruc¬ 
tion techniques are simple, but require costly equipment, and inordinate amounts 
of time. If computational reconstruction methods can be perfected, HPIV could 
be become more practical. 

There are several possibilities for computational reconstruction algorithms, 
one of which (Onural’s method) is presented in Section 5.1. The methods must 
be evaluated with realistic input in order to gauge their effectiveness. Obtaining 
realistic input is difficult. There have been experiments done using holographic 
techniques to study small ocean biological specimens, but obtaining the actual 
holograms produced is difficult. An additional difficulty is that you don’t know 
precisely what was contained in the original image so calibrating the reconstruction 
algorithm is nearly impossible. 

4.2 Hologram Synthesis 

There are several methods applicable to synthesizing holograms similar to those 
expected from HPIV. The most straightforward methods calculate the individual 
contribution from each particle to each pixel on the simulated recording plane. 
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These methods are the most accurate and flexible but suffer from excessive com¬ 
putational requirements when large numbers of particles are being simulated. 
One way around this problem is hinted at by Equation (2.18): 

= U n (x, y) <g> h Zn (; x , y) 

n 

We can subdivide the region being simulated into a large number of thin slabs. 
Within each slab there can be a large number of particles, but each will be treated 
as if it had the same 2 coordinate. We then calculate the contribution of each 
slab to the total hologram rather than the contribution from each particle. With 
a large number of particles we can achieve significant performance increases with 
little degradation of quality. 

Evaluating (2.18) is quite efficient. Using the convolution theorem: 


tp(x,y) 


^2 U n (x, y) <g> h Zn (x, y) 

n 

E ^ 1 [T[U n (x,y)\-F[KAx,y)\\ 


T~' 
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(4.1) 


Where 14, = T[U(x,y)} and 'H Zrl (ui x ,u s ) = T[h Zn (x,y)\. But H. n (a i x ,u t ) 

can be analytically evaluated [12], 


'Hzn (<*>*, w y ) = F{h Zn ) = exp 
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( w x + ul) 


(4.2) 


yielding, 


V* (a:, y) = F 


-1 


E w - (u x ,u y ) ■ exp 




(4.3) 


The algorithm to evaluate (4.3) is shown schematically in 4-1, and computer 
code to perform the algorithm is given in Sections A.1.1 and A.1.2. 
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Figure 4-1: Hologram synthesis algorithm 
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4.3 HoloSynth Simulator 



Figure 4-2: HoloSynth user interface 


4.3.1 Introduction 

HoloSynth assists in synthesizing holograms similar to those obtained from Holo¬ 
graphic Particle Image Velocimetry. It provides the user with a simple way to 
place any number of particles into a cubic volume in preparation for generating 
a synthetic hologram. The user can manually place particles anywhere in the 3- 
dimensional region as well as fill arbitrary rectangular solid regions with randomly 
distributed particles. 

HoloSynth uses the method for synthesizing holograms described in Section 
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Figure 4-3: HoloSynth main controls 

4.2. In addition it prepares an input file for an alternate hologram synthesizer 
so that different synthesis algorithms may be easily compared, without modifying 
the HoloSynth source code. The output file is also useful in analyzing the results 
of reconstruction algorithms. 

HoloSynth was developed under Red Hat®Linux version 5.1 and tested under 
Microsoft Windows NT®. 

4.3.2 The Interface 

4.3.2.1 The Main Controls 

Figure 4-3 shows the main controls for HoloSynth. 

Each button will open up a new window that controls a particular aspect of 
HoloSynth. The windows are: 

1. 3-D View of the current particle configuration. 

2. Plane View, which allows manual placement of particles on an arbitrary 
plane. 

3. Select Region, which allows selection of an arbitrary rectangular solid region 
within the view space. From this control panel you can: 

• Clear all or any part of the volume. 
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□ HoloSynth 1.0:3-D View x 



Figure 4-4: 3-D View window 


• Fill all or any part of the volume with randomly placed particles. 

4. The Synthesizer Settings window which controls the physical characteristics 
of the simulation used for the hologram synthesis. 

Through the file menu the user can open and save particle fields, and exit 
HoloSynth. 

4.3.2.2 The 3-D View 

The 3-D View, shown in Figure 4-4, provides a three dimensional view of the 
current particle configuration. There are three basic functions controlled by this 
window: the magnification (zoom) of the view, the angle and translation of the 
view, and the position of the active and two clipping planes. No direct editing 
takes place in this window, it is merely a visualization aid. 

The viewing “camera position” is controlled by the rollers and sliders to the left 
and below the viewing window. The slider above the viewing window controlling 
the magnification of the camera. 

The various view control planes are controlled using the slide bars to the right of 
the 3-D view window. They appear as colored squares drawn around the periphery 
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Figure 4-5: HoloSynth plane editor 


of the cube. Only the active plane is visible in Figure 4-4, the clipping planes are 
hidden at the ends of their travel. On a color display the active plane is bright 
green, the rear cliiping plane is blue, and the front clipping plane is magenta. The 
active plane determines the ^-coordinate that the plane editor uses when adding 
new particles. The front and rear clipping plane restrict the range of 2 coordinates 
visible in both the plane view and the 3-D view. In all cases, particles lying on 
the active plane are visible within both views. If the user wants to view only the 
particles on the active plane, then he places both clipping planes at the same limit 
of travel. 

4.3.2.3 The Plane View 

The plane view projects all particles visible in the 3-D view onto a single plane. 
Particles that are on the active (green) plane in the 3-D view are highlighted green 
in the plane view. Particles can be manually positioned or deleted on the active 
plane. 

To add a particle, simply left click within the black portion of the active plane 
view. Accurate positioning of particles can be performed in two ways: 























4.3. HOLOSYNTH SIMULATOR 


43 


n HoloSynth 1.0: 3-D Selection x 



Figure 4-6: 3-D region selection 


1. Turn on the coordinate grid. At the bottom of the plane view window is a 
selection list allowing presentation of grids of various densitites. High grid 
densities are convenient at high magnifications. While a grid is displayed 
the particles automatically snap to the nearest grid intersection. This can 
be controlled by the Display->Snap to Grid menu item, and can also be 
toggled using the ‘g’ key. 

2. Use the real-time coordinate diplay. Any time the mouse is over the edit 
window its complete 3-D coordinate is show at the bottom of the window. 

Deletion is accomplished by dragging with the right mouse button over the 
region containing the particles to be deleted. A red selection border will appear 
showing the selected region. To delete inside the selected region, use the menu 
item, Edit->Delete Inside Selection or the delete key. Using the Edit->Delete 
Outside Selection menu item will clear all particles visible in the plane editor 
but outside the selection box. The region deleted is bounded in the x and y di¬ 
rections by the selection box, and in the z dimension by the clipping planes. If 
a particle is visible in the plane editor, regardless of whether it is on the active 
plane, it can be deleted. 
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4.3.2.4 Select Regions 

The 3-D Selection window, shown in Figure 4-6, allows the user to define a solid 
rectangular sub-volume within the 3-D view and perform the following operations: 

1. Clear the entire volume. 

2. Clear all particles inside the selected sub-volume 

3. Clear all particles outside the selected sub-volume 

4. Fill the entire volume with randomly placed particles. The slider at the right 
of the window controls how many particles will be added. 

5. Fill the selected sub-volume with randomly placed particles. The slider at 
the right of the window controls how many particles will be added. 

6. Fill outside the selected sub-volume with randomly placed particles. The 
slider at the right of the window controls how many particles will be added. 

The sub-volume is displayed as a translucent yellow box with the 3-D view. Each 
wall of the selection region is controlled by its own slider. The region will only be 
displayed in the 3-D view if the Show Selection Box button has been selected. 

4.3.2.5 Hologram Synthesis 

Figure 4-7 shows the synthesizer controls. This allows easy control of the physical 
size of the simulation as well as recording plane distance and number of planes used 
to subdivide the simulation volume. The generated hologram will be displayed in a 
new window when the computations have completed. If the product is satisfactory, 
the user may save the hologram to a file using the Truevision Targa format or the 
Adobe Photshop Raw format. 

4.3.3 Requirements 

HoloSynth is written completely in C++. The Graphical User Interface uses 
the freely available platform independent user interface toolkit FLTK[13]. The 
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Figure 4-7: Synthesizer settings 


hologram synthesizer use the FFTW fast fourier transform library available freely 
from MIT[14]. The program runs under X-Windows, and Microsoft Windows with 
OpenGL[15] or the Mesa 3-D graphics library (an OpenGL functional equivalent)[16]. 

4.3.4 Program Architecture 

At its heart HoloSynth is a fancy list editor, editing a list of particle coordinates 
and radii. The core data structure is, appropriately, a linked list of particle objects. 
This allows an arbitrary number of particles to be added, restricted only by the 
memory available. 

Figure 4-8 shows the schematic relationship between all program modules that 
make up HoloSynth. The complete source code to HoloSynth is presented in 
Appendix A.l. 
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Figure 4-8: HoloSynth architecture 
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Chapter 5 

Computational Reconstruction 

5.1 Onural’s Method 

In 1985, Levent Onural proposed a technique based on a truncated inverse filter[12, 
17]. The method is summarized here, source code is presented to perform the 
algorithm in Sections A.2.1 and A.2.2. 

Starting from (2.17), 


xp z (x, y) = U (x, y) <g> h z (x, y) 

We let U (x, y) = 1 — a ( x , y) where a ( x , y) is the opacity of the given plane. Now 
write the recording plane intensity, 

I z (x,y) = i/j z (x,y)ip* z (x,y) 

= |(1 a (a;, t/)) <S> h z (a;, y)| 

= 1 - a* (x, y) <g> h* z (x, y) - a (x, ?/) <g> h 2 (x, y) 

+ \a(x,y) ®h z (x,y)\ 2 
= 1 -a(x,y)® [^(x,?/) + /i 2 (x,y)] 

= 1 - a(x,y) ® 2Re{h x (x,y)} 

Where \a(x,y) ®h z (x, y)\ 2 was dropped as insignificant in (5.1). Equation (5.2) 


(5.1) 

(5.2) 
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shows the intensity modeled as a DC shifted linear system with impulse response, 

9 fay) = 2Re{h z (x,y)} (5.3) 

No inverse of g(x,y) in the normal sense exists, However, Onural showed in 
[12] that a series approximation to the Fourier transform of g can be found. If 
Q (u x ,ujy) = T[g (x, y)\ then, 


r>—\ t Xz ( 2 , 2\ , ^ f M 1 Jlo§2 1 o 

Gm fax, Uy) - COS — (u x + U y ) + 2^ < -TV- f S, 

L—1 V ' 


kfax,^y) (5-4) 


where [log 2 k} is the integer obtained by truncating the value of log 2 k and, 


S k (u x , uj y ) = (-l) fc cos ^ (2k + 1) (ul + ujy) 


yielding, 


9M fa’V) = ? 1 [Gm fax,U y j\ 


Looking back at the hologram recording given by (5.2), 

a ( x , y) <g> 2 Re {h z (x, y)} = l- I (x, y) 
T[a(x,y)\g(u x ,u y ) = F[l - I (x,y)] 


a z(x,y ) — T [^[1 I (x,y)]QjJ z (u x ,Uy)^ (5.7) 


yielding the original opacity function. 


5.2 Simulation Results 

5.2.1 Synthesis and Computational Reconstruction 

Figure 5-2 is a portion of a hologram generated by HoloSynth. The particle config¬ 
uration consisted of 2328 particles. 2000 of the particles were randomly distributed 
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Figure 5-1: Schematic of the particle configuration used to generate Figure 5-2 

throughout a 5cm cube. Along the center plane parallel to the wavefront 328 par¬ 
ticles were positioned to spell the letters “MIT”. Figure 5-1 illustrates a schematic 
of that configuration. The synthetic hologram itself was 2048x2048 pixels, and 
the volume was subdivided into 1024 planes, for a 2-axis discretization of 0.04mm. 
The simulated particle diameter was 25/rm. The portion of the hologram shown in 
5-2 is 512x512 pixels and its upper left corner is 205 pixels from the left extremity 
of the original and 741 pixels from the top. Although it contains the C M’ from 
“MIT” , it is impossible to distinguish from the diffraction patterns produced by 
the other particles. 

Figure 5-3 through 5-8 show the results of reconstructing Figure 5-2 at various 
distances between 30.0cm and 30.5cm from the recording plane using OnuraPs 
method with M — 5. The images generated are the opacity functions, with black 
being transparent, white being opaque. 

In Figure 5-3(a) the ‘M’ is clearly visible on a computer monitor with some 
slight clouds from out-of-focus particles. Figure 5-4 is the same hologram recon¬ 
structed at 30.1cm. Here there is little visible difference from the previous figure. 
This is an artifact of the way the decoder transform real valued intensities into 
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gray scale bitmap images. 

To translate from an arbitrary field of real number to a field of 8 bit integers, 
the decoder first finds the maximum and minimum values present in the real array, 
it then scales and translates each of these values to fill the range from 0 to 256, 
in order to maximize use of dynamic range available in an 8 bit gray scale image. 

The 30.1cm reconstruction has a much lower maximum intensity than the 
30.0 cm reconstruction. Volume reconstruction, where many planes will be recon¬ 
structed at a time, in order to reconstruct flow fields, will have to keep track of 
this scaling issue between planes. 

The important thing to note about this series of decodes, is that while the 
source hologram was constructed using a 2 discretization of 0.04mm, the recon¬ 
struction algorithm 2 -discrimination was of the order of millimeters. Thus, the 
original image could have been discretized into 64 steps, or 0.78mm slabs, requiring 
1 /16th the computation time with little or no loss of reconstruction accuracy. 

5.2.2 Image Postprocessing 

The next step in reconstruction is to automatically detect particles in each plane 
and register their x, y, and 2 coordinates, rather than simply displaying an im¬ 
age. Repeatedly reconstructing at different depths and detecting the particles will 
replace the time consuming 3-D scanning. 

In Figures 5-3 through 5-8 subfigure (b) and (c) show the results of applying 
some simple image processing techniques to the original reconstructions to enhance 
the output. Subfigure (b) in each figure applies ‘edge detection’, which simply 
filters the image marking any region where the gradient of the image exceeds 
a certain threshhold. In the case of the grayscale images being used here, the 
gradient is simply the difference in grayscale value between adjacent pixels. [18] 
Subfigure (c) in each figure ‘threshholds’ the image, changing a pixel to white 
(value 255) or black (value 0) depending on whether or not the pixel value exceeds 
a preset threshhold. 1 


’ The image processing was accomplished with the freely available GNU Image Manipulation 
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The combination of edge detection and threshholding can be used very effec¬ 
tively to automatically detect particle positions. At the end of the image process¬ 
ing, there are only regions of black (value 0) with very small regions of white (value 
255). The white regions correspond to particles and can easily be automatically 
tabulated and mapped into x and y coordinates. 


Program (GIMP, http://www.gimp.org), which is very similar to Adobe Photoshop®. 
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Figure 5-2: Portion of a synthetic hologram. The original synthetic hologram was 
2048x2048 pixels. This portion is 512x512 pixels, cropped from the orig¬ 
inal hologram to make the images presentable in printed form. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de¬ 

tection. 



(c) Image after enhancing with edge de¬ 
tection and threshholding. 


Figure 5-3: Figure 5-2 decoded at a depth of 0.300m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with' edge de¬ 

tection. 



(c) Image after enhancing with edge de¬ 
tection and threshholding. 


Figure 5-4: Figure 5-2 decoded at a depth of 0.301m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de¬ 

tection. 



(c) Image after enhancing with edge de¬ 
tection and threshholding. 


Figure 5-5: Figure 5-2 decoded at a depth of 0.302m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de¬ 

tection. 



(c) Image after enhancing with edge de¬ 
tection and threshholding. 

Figure 5-6: Figure 5-2 decoded at a depth of 0.303m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de¬ 

tection. 



(c) Image after enhancing with edge de¬ 
tection and threshholding. 


Figure 5-7: Figure 5-2 decoded at a depth of 0.304m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de¬ 

tection. 



(c) Image after enhancing with edge de¬ 
tection and threshholding. 


Figure 5-8: Figure 5-2 decoded at a depth of 0.305m. 
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Chapter 6 

Conclusions and 
Recommendations 


It is clear that computational reconstruction of HPIV is practicable. The proce¬ 
dure that Onural proposed in 1985 was not then practical due to the computational 
power commonly available. Modern computing can easily handle realistic holo¬ 
grams and the numeric techniques of normal PIV could be extended to deal with 
the particle identification problem in HPIV. The simulator presented in this thesis 
will be useful in any setting when attempting to validate particular reconstruction 
algorithms. 

Future work should concentrate on extending the reconstruction prosses fully 
into three dimensions and developing particle identification techniques. Automat¬ 
ically determining the velocity vector for a particle is the ultimate goal. 
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Appendix A 


Source code 


A.l HoloSynth: Hologram Synthesizer. 

This is the source code for HoloSynth, a program for easily creating particle configuration 
and holograms for use calibrating computational reconstruction algorithms. 

A. 1.1 OnuralSynthesizer.hpp 

#ifndef ONURALSYNTHESIZER_H 

// {include guard 

#def ine ONURALSYNTHESIZERJH 1 

#include ’’Synthesizer.hpp” 

#include ’’ParticleList .hpp” 

#ifndef FLMATH 
#define FLMATH 
#include <FL/math.h> 

#endif 

#include ’’RGBImage .hpp” 

/ * !\ class OnuralSynthesizer 
* 

*\brief The algorithm use for actually generating the hologram 
* 

* The forward diffraction algorithm implement a method described in 

* ’’Digital Decoding of In-Line Holograms” a PhD thesis by Leven Onural 

* */ 

class OnuralSynthesizeripublic Synthesizer 

{ 

public: 

/// created with a reference to the current particle list 
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OnuralSynthesizer(ParticleList *pl); 

~ OnuralSy nt hesizer (); 

/// Generate, organized all the FFTs and diffraction calculations. 
int generate(); 

private: 

RGBImage *hologram; 

void diffract (fftw.complex *image); 
void clearField(fftw_complex *field, fftw_real val); 
void plot Particle (particle *c,fftw.complex *plane); 
void plot (int pJt, int p_Y, fftw_complex *plane); 
void addto(fftw.complex *al, fftw.complex *a2); 
void MagtoRGB(fftw_complex *data); 

void cleanup(); 

fftwnd.plan fftwFwdPlan; 
fftwnd.plan fftwBwdPlan; 

fftw.complex *oplane; 
fftw.complex *hplane; 


#endif 

A.1.2 OnuralSynthesizer.hpp 

#include <FL/F1.H> 

#include’’ProgressUI.hpp” 

#include ’’ScrollImageUI .hpp” 

#include <FL/Fl_Image. H> 

#include ’’OnuralSynthesizer.hpp” 

#include ’’StatusBox.hpp” 

#define MIN(a, b) (((a) < (b)) ? (a) : (b)) 
#define MAX(a, b) (((a) > (b)) ? (a) : (b)) 


OnuralSynthesizer::OnuralSynthesizer(ParticleList *pl) 
:Synthesizer(pl) 

{ 

} 

OnuralSynthesizer::~OnuralSynthesizer() 

{ 


cleanupQ; 
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} 

int OnuralSynthesizer::generate() 

{ 

int progresscount=0; 

progressUI—>-setMax(5+5*dJPlanes); 
progressUI—Hipdate(O); 
progressUI^reset(); 
progressUI ->show (); 

Fl::wait(0.0); 

particle *curr; 

int pixels, il, pcount, Neale, Mcalc; 
fftw_real zgrain, Zmin, Zmax; 


* Initialization. 

pixels==2*(2*N/2+l)*2*M; 

Ncalc=2*N; / /Double the size for padding out the edge effects. 
Mcalc=2*M; 

if(progressUI—!>setStatus( "Allocating memory. . . ")){ 
cleanup(); 
progressUI—>*hide(); 

return 0; 

}; 

hplane=(fftw_complex *) calloc( pixels,sizeof(fftw_complex)); 
oplane=(fftw_complex *) malloc( pixels*sizeof(fftw_complex)); 

progressUI—Mipdate(++progresscount); 

zgrain=(fftw_real)l/(2*dJPlanes);//77ie thickness of a HALF plane 

if(progressUI-*setStatus(" Generating FFTV Plans...")){ 
cleanup(); 
progressUI—^hide(); 

return 0; 

}; 


fftwFwdPlan=fftw2d_create_plan(Ncalc, Mcalc, FFTW_FORWARD, 
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FFTWJN.PL ACE); 

fftwBwdPlan=fftw2d_create_plan(Ncalc, Mcalc, FFTW_BACKWARD, 
FFT W JN .PLACE); 

progressUI->update(++progresscount); 

//Remember we will eventually need to split out the vector end of things. 
pList—>sortList OnZ (); 


curr=pList—»first(); 
for(il=l;il<=d_Planes;il++){ 


Fl::wait(0); 

if(progressUI—>setStatus(" Clearing Object Plane... ")){ 
cleanup(); 
progressUI—>hide(); 

return 0; 

}; 

clearField(oplane, 1.0); 
progressUI*->update(++progresscount); 

Zmin=(il-l)*2*zgrain-0.5; 
Z=(Zmin+zgrain)*depth;//P%S 2 caZ Z 
Zmax=Zmin+2*zgrain; 

pcount=0; 

if(progressUI—)>set Status ("Plot ting Particles. . .")){ 
cleanup (); 
progressUI-->hide(); 

return 0; 

}; 


while((curr!=NULL)&;&;(curr—►Z>Zmin)&;&;(curr—>Z<=Zmax)){ 
plotParticle(curr,oplane); 
pcount++; 
curr=curr—>next; 

} 

progressUI—>-update(++progresscoimt); 
if(pcount>0){ 


if(progressUI—»set,Status( "Performing Forward FFT...")){ 



A.l. HOLOSYNTH: HOLOGRAM SYNTHESIZER. 


64 


cleanup(); 
progress!!!—»hide(); 

return 0; 

}; 

fftwnd_one(fftwFwdPlan, oplane, 0); 
progressUI—>update(++progresscount); 

if(progressUI—^setStatus("Diffracting... ")){ 
cleanup(); 
pr ogressUI-*hide (); 

return 0; 

}; 

diffract (oplane); 

progressUI—mpdate(++progresscount); 

if(progressUI-^setStatus(" Adding...")){ 
cleanup(); 
progressUI-^hideQ; 

return 0; 

}; 

addto (hplane, oplane); 

pr ogr essUI —> up dat e (■++pr ogr esscount); 

} 

progresscount=5*il+2; 
progressUI—mpdate(progresscount); 

} 

hplane [0] .re= 1000000; 

progressUI—»setStatus("Performing Backward FFT... n ); 
progressUI—>update(++progresscount); 
fftwnd_one(fftwBwdPlan, hplane, 0); 
progressUI—)'SetStatus( ,, Calculating Intensity. .. "); 

MagtoRGB (hplane); 

//The hologram is finished, now display it.. 

FIJmage *image=new FlJmage(hologram“»data(), N, M); 
progressUI—»setStatus( "Displaying.. . "); 
progressUI—mpdate(++progresscount); 

ScrollImageUI *si=new ScrollImageUI(image , 0,0, 500,500); 
si^show(); 

progressUI-^hide(); 

cleanup(); 


} 


return 1; 
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void OnuralSynthesizer::addto(fftw_complex *a, fftw.complex *b) 

{ 

int il; 

for (il=0; il < ( 4 * N * M ); il++){ 

if((fabs(b[il].re)>=0.0000001)&&(fabs(b[il].im)>=0.0000001)) 

{ 

a[il].re += b[il].re; 
a[il].im += b[il].im; 

} 

} 

} 

void OnuralSynthesizer::clearField(fftw_complex *field, fftw_real val) 

{ 


int il; 

for(il=0; il < ( 4 * N * M ); il++){ 
field[il].re = val; 
field[il].im = 0.0; 


} 

void OnuraISynthesizer::plotParticle(particle *c,fftw_complex *plane) 

{ 

int x, y, r; 


x=(int)((c—+ 0.5)*N); 
y=(int)((0.5 - c-*Y)*M); 
r=(int)(c—>-R*N); 


switch (r) 

{ 

case 6: 

plot(x-l, y+1, plane); 
plot(x-l, y-1, plane); 
plot(x+l, y-1, plane); 
plot(x+l, y+1, plane); 
case 2: 
case 3: 
case 4: 
case 5: 


plot(x, y-1, plane); 
plot(x, y+1, plane); 
plot(x-l, y, plane); 
plot(x+l, y, plane); 
default: 

plot(x, y, plane); 
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} 


} 

void OnuralSynthesizer::plot(int p_X, int p_Y, fftw_complex *plane) 

{ 

int pos=2 * p_Y * N + p_X; 
if(pos>0){ 

plane[ pos].re=0.0; 
plane[ pos].im=0.0; 

} 


} 

#def ine FJPI 4.0*3.1415926535 

void OnuralSynthesizer::diffract(fftw_complex *image) 

{ 

int x, y, Neale, Mcalc; 
long ik; 
fftw_real lz, cc; 
fftw-real wx, wy, a, ca, sa; 

Ncalc=2*N; 

Mcalc=2*M; 

lz=L* (distance-Z); 
ik=0; 

/* FFTW puts frequency space origin at the upper right corner of the 

* array. Positive frequencies build from the upper left comer 

* diagonally down through the upper left quadrant The other corners 

* of the array the array are ALSO at the origin, so the very center 

* of the array is the HIGHEST frequency. 

* 

* In order to use the diffraction filters we MUST take that into 

* account, that is what the strange trigraph statements take care of 

* below */ 

for (y=0;y<Mcalc;y++) 

{ 

wy = (y<M) ? y/height : (y-Mcalc)/height; 

// wy — (y<M) ? y/height : (Mcalc-y)/height; 
wy * = wy; 

for(x=0;x<Ncalc;x++){ 

wx=(x<N) ? x/width : (x-Ncalc)/width; 

// wx=(x<N) ? x/width : (Ncalc-x)/width; 

a=lz*(wy+wx*wx)*3.1415926535; 

// a=lz*(wy+wx*wx)/F~PI; 
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ca=cos(a); 

sa=sin(a); 

cc=image[ik] .re*ca-image[ik] .im*sa; 
image[ik].im=image[ik].re*sa+image[ik].im*ca; 
image[ik].re=cc; 
ikH—f; 

} 

} 

} 

void OnuralSynthesizer::MagtoRGB(fftw_complex *data) 

{ 

int i, j, ij, grey, Neale, Mcalc; 
fftw_real v, vmax, vmin, vl, v2; 


Ncalc=2*N; 

Mcalc=2*M; 

hologram ==new RGBImage(N,M); 

vmax=0.0; 

vmin=le6; 

ij=0; 

/ / Scale the result to take maximum, advantage of 8bit range 
for(i=0;i<Mcalc*Ncalc;i++){ 
vl=data[i].re; 
v2=data[i].im; 
v=sqrt(vl*vl+v2*v2); 

vmax=MAX(v, vmax); 
vmin=MIN(v, vmin); 

} 

vmax=l / (vmax-vmin); 


RGB Color c; 
for(i=0;i<M;i++){ 
for(j=0;j<N;j++){ 
ij=i*Ncalc+j; 


vl=data[ij].im; 

v2=data[ij].re; 

v=(sqrt (■vl * vl+v2* v2)-vmin) * vmax; 
data[ij].re=v; 
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grey=(int)floor((fftw_real)(v*255.0)); 

c.B=grey; 

c.R=grey; 

c.G=grey; 

hologram—>plot(j, i, c); 

} 

} 


} 

void OnuralSynthesizer: :cleanup () 

{ 

free(hplane); 

free(oplane); 

fft,wnd_destroy_plan(fftwF wdPlan); 
fftwnd_destroy_plan (fftwBwdPlan); 

} 
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A.2 Decoder: Onural Reconstruction Algorithm. 


This is the source code for an Onural reconstruction algorithm. 


A.2.1 OnuralDecoder.hpp 


#ifndef ONURALDECODERJiPP 
#def ine ONURALDECODER.HPP 1 
#include ’’Decoder .hpp” 

#include ’’ProgressUI.hpp” 

#include <fftw.h> 

class OnuralDecoder: public Decoder 


{ 


public: 

OnuralDecoder (); 

virtual uchar* decode(uchar *data, int N, int M, 

float W, float H, float L, 
float D); 


private: 

int dJmageN; 
int d-CalcN; 
int dJmageM; 
int d.CalcM; 

float dJmageW; 
float dJmageH; 
float dJmageLambda; 
float d_decodeDist; 

fftw.complex *d_Field; 

uchar * d_data; 


void OnuralDecoder::noDClquad(fftw-complex *field, fftw_real average Value); 

fft,w_real OnuralDecoder: :aprinv(fftw_real zz, int Mstage); 

void OnuralDecoder::onural(fftw-complex *hg, 

int N, int M, 

fftwjreal L, fftw_real Z, 

fftwjreal X, fftw_real Y,int Mstage); 


fftw-real Magnitude(fFtw_complex *field); 

void uchar2field(); 
uchar *field2uchar(); 
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fftw-real *clipField(); 


uchar *OnuralDecoder: :mainLoop (); 
ProgressUI ^progress; 


}; 

#endif 


A.2.2 OnuralDecoder.cpp 

#include ’’OnuralDecoder .hpp” 

#include <fftw.h> 

#include ’’StatusBox .hpp” 

#include <FL/F1.H> 

#include <math.h> 

#define PI 3.14159265359 
#def ine TPI 6.28318530718 
#define MIN(a, b) (((a) < (b)) ? (a) : (b)) 
#define MAX(a, b) (((a) > (b)) ? (a) : (b)) 

OnuralDecoder::OnuralDecoder() 

{ 

dJmageN=0; 

dJmageM=0; 

dJmageW=0; 

dJmageH=0; 

dJmageLambda=0; 

d_decodeDist=0; 

d_data=0; 

progress =new ProgressUI(); 


uchar *OnuralDecoder::decode(uchar* data, int N, int M, float W, float H, float L, 
float D) 

{ 

dJmageN=N; 

d_CalcN=2*N; 

dJmageM=M; 

d_CalcM=2*M; 

dJmageW=W; 

dJmageH=H; 

dJmageLambda=L; 

d_decodeDist=D; 

d_data=data; 
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return mainLoopQ; 


/* Remove the DC component from only the first Quadrant, the other 
* three quadrants are ALREADY ZERO*I 

void OnuralDecoder::noDClquad(fftw_complex * field. fftw_real average Value) 

{ 

int i, j, pos; 

for(i=0; i<d_ImageM; i++){ 

for(j=0; j<dJmageN; j++){ 
pos=i*d_CalcN+j; 
field[pos].re -= averageValue; 

} 

} 


fftw_real OnuralDecoder::aprinv(fftw_real zz, int Mstage) 

{ 

fftw_real cc,ss,powKa; 
int k,powK; 

cc=cos(zz); 

ss=0; 

if(Mstage>l){ 
ss=2*sin(zz); 
cc+=ss*sin(2*zz); 
ss* =2*cos(2*zz); 

} 

for (k=3;k<=Mstage;k++){ 

powK=4<C(k-3); // powK=2**(k-l), fine for small k, and FAST 

powKa=(fftw_real)powK*zz; 

cc+=ss*sin(powKa); 
ss* =2*cos(powKa); 

} 

return 2*cc/Mstage; 

} 

void OnuralDecoder::onural(fftw_complex *hg, 
int N, int M, 
fftw_real L, fftw_real Z, 
fftw_real X, fftw_real Y,int Mstage) 
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int x, y, ml2, nl2; 

long ik; 

fftw_real alpha, lz, wx2, wy2; 

lz=L*Z; 

ml2=M/2; 

nl2=N/2; 

ik=0; 


for (y=0;y<M;y ++)//M 

{ 

wy2 = (y<ml2) ? y/Y : (y-M)/Y; 
wy2 * = wy2; 
for(x=0;x<N;x++)//7V 
{ 

wx2 = (x<nl2) ? x/X : (x-N)/X; 
wx2 * = wx2; 

alpha=PI*lz* (wy2+wx2); 

hg[ik].re* =aprinv(alpha, Mstage); 
ikH—h; 

} 

if(progress-^update(3)){ 
progress-* hide(); 

return; 

}; 

} - 

} 

fftw_real OnuralDecoder::Magnitude(fftw_complex * field) 

{ 

int i; 

fftw_real v,sumV, vl, v2; 
sumV=0; 

for(i=0; i < d_CalcN*d_CalcM; i++){ 
vl=field[i].re; 
v2=field[i].im; 
v=sqrt(vl*vl+v2*v2); 
field [i].re=v; 
field[i].im=0; 
sumV+=v; 
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} 


return sumV/(cLCalcN * d.CalcM); 


} 

void OnuralDecoder::uchar2field() 

{ 

int arraysize=d_CalcN*d_CalcM*sizeof(fftw_complex); 
int fieldpos, ucharpos, x, y; 
uchar minval, maxval; 
fft,w_real scale; 

dJField=(fFtw_complex *)calloc(arraysize, 1 );//calloc clears to zero ... 

minval=0; 

maxval=255; 

for(y=0; y< dJmageM; y++ ){ 
for(x=0; x<dJmageN; x++){ 

ucharpos=3*(x -j- y*dJmageM); 
minval=MIN(minval, d_data[ucharpos]); 
maxval=MAX(maxval, d_data[ucharpos]); 

} 

} 

scale= 1.0/ (fftw_real) (maxval-minval); 

for( y=0; y< dJmageM; y++ ){ 
for( x=0; x<dJmageN; x++){ 

fieldpos=x + y*d_CalcN; 

_ucharpos=3*(x + y*dJmageN); 

d-Field [fieldpos] .re= ((fftw_real) (d_data[ucharpos]-minval)) *scale; 

} 

} 


} 

uchar *OnuralDecoder::field2uchar() 

{ 

int arraysize=3*dJmageN*d_ImageM*sizeof(uchar); 
int fieldpos, ucharpos, x, y; 
fftw_real minval, maxval, v, vre, vim; 
fFtw_real scale, val; 

uchar *retval; 

retval= (uchar* )calloc(arraysize, 1 )\//calloc clears to zero... 
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minval==le20; 

maxval=0; 

for(y=0; y< dJmageM; y++ ){ 
for(x=0; x<dJmageN; x++){ 
fieldpos=x + y*d_CalcM; 
vre=dJField[fieldpos].re; 
vim=dJField[fieldpos].iin; 
v=sqrt (vre* vre+vim* vim); 

minval=MIN(minval, v); 
maxval=MAX(maxval, v); 
d -Field [fieldpos]. re=v; 


} 

} 

scale=255.0/(maxval-minval); 

for(y=0; y< dJmageM; y++ ){ 
for(x=0; x<dJmageN; x++){ 
fieldpos=x + y*d_CalcN; 
ucharpos=3*(x + y*dJmageN); 

val=(uchar)((d-Field[fieldpos].re-minval)*scale); 
ret val [uchar p os]—(uchar) val; 
retval[ucharpos+l]=(uchar)val; 
ret val [uchar p os+2]=(uchar) val; 


} 

} 

return retval; 


uchar* OnuralDecoder::mainLoop() 


fftw_real DC; 

int progresscount=0; 

progress-»setMax(5); 
progress—mpdate(O); 
progress—> reset (); 
progress—>show(); 

if(progress-»setStatus( M Transforming Input. ..")){ 
progress—>hide(); 
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return 0; 

}; 


fftwnd.plan ForwardFFTPlan, BackwardFFTPlan; 
uchar2field(); 

progress—mpdate(+H-progresscount); 
if(progress—^setStatus(" Making FFTW Plans... ")){ 
progress—>hide(); 

return 0; 

}; 

ForwardFFTPlan=fftw2d_create_plan(d_CalcN, d.CalcM, 

FFT W.FORWARD, 

FFTW .ESTIMATE | FFTW JN .PLACE); 

BackwardFFTPlan=fftw2d_create_plan(d_CalcN, d.CalcM, 

FFTW-BACKWARD, 

FFTW .ESTIMATE |FFTW JN .PLACE); 

DC=4*Magnitude(d_Field); //The 4 handles the doubling of array size 


noDClquad(d_Field, DC); 

progress—»update(++progresscount); 

if(progress—»setStatus( "Calculating Forward transform..."))! 
progress—>hide(); 

return 0; 

}; 

fftwnd_one(ForwardFFTPlan, d_Field, 0); 

d_Field[0].re=0; 

d_Field[0].im=0; 

progress—>update(++progresscount); 
if(progress—>setStatus(" Starting Onural reconstruction..."))! 
progress—»hide (); 

return 0; 

}; 

onural(d_Field, d.CalcN, d.CalcM, dJmageLambda, d_decodeDist, 
dJmageW, dJmageH, 5); 
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progress —> up dat. e (■+ •4- pr ogr esscount); 
if(progress—»setStatus( "Calculating the Reverse Transf orm. . . ")){ 
progress—>Lide(); 

return 0; 


fftwnd_one(BackwardFFTPlan, d_Field, 0); 

progress-»update(++progresscount); 
if(progress-*setStatus(" Complete ")){ 
progress—>*hide(); 

return 0; 


uchar* decoded=field2uchar(); 

free(dJField); 

progress—^hide(); 
return decoded; 


}; 
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